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Level-set topology optimization is used to design a wing considering skin buckling under 
static aeroelastic trim loading, as well as dynamic aeroelastic stability (flutter). The level-set 
function is defined over the entire 3D volume of a transport aircraft wing box. Therefore, the 
approach is not limited by any predeflned structure and can explore novel conflgurations. 
The Sequential Linear Programming (SLP) level-set method is used to solve the constrained 
optimization problems. The proposed method is demonstrated using three problems with 
mass, linear buckling and flutter objective and/or constraints. A constraint aggregation 
method is used to handle multiple buckling constraints in the wing skins. A continuous 
flutter constraint formulation is used to handle difflculties arising from discontinuities in the 
design space caused by a switching of the critical flutter mode. 


I. Introduction 

T opology optimization is the most general form of structural optimization, as it is the least restricted by a pre- 
determined design [1]. Therefore, topology optimization is an ideal tool for exploring novel structural 
configurations that may be far from the conventional designs. For aircraft wing design, an efficient structure is one 
that operates at minimum mass, whilst meeting a series of constraints, which may include skin buckling, stress and 
dynamic properties, such as natural frequencies and flutter. Topology optimization has been applied to improve 
aircraft wing structure designs. Some approaches only consider a simple volume constraint, whilst optimizing 
objectives such as compliance [2-5], aileron reversal speed [6], natural frequencies [7] and aeroelastic stability [7, 
8]. The work of Eves et al. [9] also included tip twist and extrusion manufacturing constraints with the usual volume 
constraint in a minimization of compliance problem. Optimization for skin buckling and stresses were then 
considered in separate steps from the topology optimization of the wing structure, which was mainly used to guide 
the placement of spar-like members. Eigenvalue constraints have also been used when optimizing for the aeroelastic 
[10] or aerothermoelastic [11] performance of wing structures. James et al. [12] used topology optimization in a 
multi-disciplinary framework to optimize the aerodynamic performance of a wing under static aeroelastic loading. 
The objective was to minimize drag whilst meeting several lift constraints. Maute and Allen [13] minimized wing 
mass, subject to structural and aerodynamic constraints, including tip displacement, skin stresses, wing stiffness, lift 
and drag. This approach also considered the aeroelastic coupling between the wing deflection and static 
aerodynamic loading. The aeroelastic coupling of static aerodynamic loading and deformed wing shape is an 
important consideration in the topology optimization of aircraft wings, as the coupling effect can strongly influence 
the optimal design [12, 13]. 
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There are several approaches to structural topology optimization. The most popular approach applied to aircraft 
wing design is the Solid Isotropic Material with Penalization (SIMP) method [I]. However, this approach can 
produce optimal designs with a significant amount of intermediate density material, which can be difficult to 
interpret as a manufacturable design. This is particularly noticeable when the problem includes multiple constraints 
and the aeroelastic coupling effect is taken into account during optimization [12, 13]. An alternative is to use the 
level-set topology optimization method, which defines the structure using an implicit function and produces 
solutions with clear boundaries with no intermediate density material. The implicit description of the structure 
naturally allows for topological changes, such as merging or splitting of holes. Several level-set based optimization 
methods have been developed in recent years [14-16]. The conventional approach is to define the structure using an 
implicit signed distance function, which is updated by solving a Hamilton- Jacobi (H-J) type equation, where the 
velocity function is derived from the shape derivatives of the objective and constraints. 

A key feature of aeroelastic wing optimization is a potentially large number of design constraints and heavily 
constrained problems have historically proved to be challenging for level-set methods. A number of approaches 
have been proposed to perform constrained optimization using the conventional level-set method, which are 
summarized here. Early methods used a Lagrange multiplier to add a single constraint to the objective, thus turning 
the constrained problem into an unconstrained one. Different methods have been used to compute the Lagrange 
multiplier, including the fixed value approach [17] (equivalent to a penalty method) and using a one dimensional 
line search to find a value that maintains constraint feasibility, or moves the design towards the feasible region [18]. 
Lor problems with multiple constraints a gradient projection method has been used [19], although a return-mapping 
algorithm is required to handle infeasible designs and non-linear constraints. The augmented Lagrange multiplier 
method has also been used to handle multiple constraints [20], although the choice of penalty parameters can affect 
the efficiency and convergence behavior of the level-set method [21]. A recent development is the Sequential Linear 
Programming (SLP) level-set method [22]. This approach uses discretized boundary integrals of the shape 
derivatives and a SLP sub-problem to compute a velocity function that improves the objective, whilst meeting the 
constraints. This method can also simultaneously optimize for additional design variables not explicitly related to 
the level-set method, such as sizing parameters and support boundary conditions. 

There have been few applications of the level-set method for wing optimization. Gomes et al. [6] used the 
spectral level-set method to optimize skin reinforcement for enhanced roll maneuvers. James et al. [4] used an 
isoparametric level-set method to minimize the compliance of a 3D wing box subject to static aerodynamic loads. In 
our previous work the level-set method was developed to minimize the compliance of a realistic 3D wing box under 
static aeroelastic and self-weight body force loading [23]. A 3D unstructured level-set method was developed to 
optimize the entire volume of the 3D wing box and thus the solution was not restricted by a predefined structural 
layout. Other applications of topology optimization for aircraft wing structure design are often limited to local 
regions or single components, such as ribs [3, 24], stiffeners [11] or skins [6], or it is used to improve an existing 
structural layout that resembles the classic spar-rib configuration [7, 13]. The SIMP method has also been applied to 
optimize the full wing box volume. However, the approaches by Eves et al. [9] and Oktay et al. [5] did not include 
the aeroelastic coupling and the solutions obtained by James et al. [12] contained significant intermediate density 
material. 

In our previous work [23], the 3D wing box design space was discretized using 3D solid finite elements. 
However, the wing box must have at least an upper and lower skin to support the aerodynamic loads, and will also 
require leading and trailing edge spars, to provide attachment points for leading and trailing edge devices, to form a 
box to contain the fuel and to provide some connectivity between the upper and lower skins. These features can be 
modeled by fixing 3D solid finite elements along the boundary of the wing box to remain part of the structure. 
However, if constraints are applied to these regions, for example skin buckling, then the solid element modeling of 
the thin shell type regions can lead to a poor approximation of the constraint function. An alternative approach is 
used in [9], where an outer layer of shell elements are connected to the solid element mesh used for topology 
optimization. In this paper we adopt a similar strategy where the upper and lower skins and leading and trailing edge 
spars are modeled using shell elements. The thickness of these elements can be fixed or treated as additional design 
variables in the optimization problem by using the SLP level-set method. Including shell thickness design variables 
expands the design space and therefore increases the potential for finding an improved design. 

The goal of this paper is to use the SLP level-set topology optimization method to design an aircraft wing 
considering skin buckling under static aeroelastic loads and dynamic aeroelastic stability (i.e. flutter). Three 
problems are considered with mass, flutter dynamic pressure and linear buckling as the objective or constraints. The 
effect of the static aeroelastic loading on the buckling load factor is taken into account through a multi-disciplinary 
analysis of the system and a coupled adjoint sensitivity analysis, flutter sensitivities are computed with derivative 
equations for nonlinear eigenvalue migration problems, and multiple instability failure modes are aggregated into a 
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single continuous formulation. Methods are presented to calculate the shape derivative of the aeroelastic objectives 
and constraints, which are required by the SLP level-set method to obtain the velocity function used to optimize the 
level-set function. 


II. Level-Set Topology Optimization 

This section briefly outlines the SLP level-set method and the schemes used to perform level-set optimization on 
3D unstructured grids. The SLP level-set method is based on the conventional level-set method. First, the boundary 
of the structure is defined as the zero level-set of an implicit function: 


(p[x^ > 0,x G Q 

^ ^^(x) = 0,xgL (1) 

^(x) < 0,x ^ Q 


where ^x) is the implicit function and x g where Q.d is the design domain containing the structure, Q The 
implicit description of the structure allows the boundary to break and merge naturally, enabling topological changes 
to occur during optimization. 

The level-set description of the structural boundary can be used for optimization. First, a velocity field is derived 
from shape sensitivity analysis. The implicit function is then updated by solving the equation: 


dt 


+ V^(x,/)— = 0 


( 2 ) 


where Ms a fictitious time domain. Equation (2) becomes a Hamilton-Jacobi type equation by considering advection 
velocities of the type: dx ! dt = V • (x). The equation is then discretized and solved using an explicit forward 

Euler scheme, producing the following update rule for the discrete implicit function values: 

(3) 


where i is a discrete point within the domain, k is the current iteration number, E is a velocity function defined 
normal to the boundary, such that a positive velocity moves the boundary inwards, and is a discrete time step. 

To perform optimization using the level-set method, the velocity function in Eq. (3) is often defined from shape 
derivatives. Shape derivatives usually take the form of a boundary integral of a shape sensitivity function multiplied 
by the velocity function: 



( 4 ) 


where / is a function dependent on the structural domain and Sf is the shape sensitivity function for /, which is 
continuous along the boundary. A velocity function is sought that changes the structural domain such that the 
objective is improved and the constraints remain feasible, or move towards the feasible region. In conventional 
level-set optimization the design variable is effectively the structural domain Q, which is implicitly defined by the 
level-set function, Eq. (1). 

The SLP level-set method is used to obtain the velocity function at the boundary interface to enable constrained 
optimization. Full details of the method can be found in [22]. A key aspect of the SLP level-set method is the 
discretization of the shape derivative boundary integrals (e.g. Eq. (4)) for the objective and constraint functions. This 
provides an estimate of how each function will change when the implicit function is updated under a velocity 
function. The velocity function is then obtained by solving a sub-problem using SLP. The sub-problem is formulated 
to maximize the improvement of the objective function, whilst maintaining the feasibility of the constraint functions, 
or at least moving the design towards the feasible region. The step change for non-level-set variables, such as wing 
skin thickness, can be added to the SLP sub-problem, enabling simultaneous optimization of topology and sizing 
variables. 
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We now briefly describe the implementation of the conventional level-set method for unstructured grids. For more 
details see [23]. The implicit hmction is discretized on a grid that is fitted to the domain. For convenience this grid is 
also used for the finite element analysis (FEA) although this is not necessary. To avoid expensive re-meshing, the 
FEA mesh remains fixed throughout the optimization. This produces elements cut by the structural boundary, as 
defined by the implicit function Eq. (1). The properties of a cut element are approximated using the volume-fraction 
weighted method, where stiffness and mass are scaled by the fraction of the element that lies within the structure. 
Discrete values of the shape sensitivity function in Eq. (4) are computed at the boundary interface using a weighted- 
least squares method [25], which reduces the errors introduced by the volume-fraction approximation at the 
boundary. The velocity function obtained using the SEP level-set method is only specified at the boundary interface. 
To update the implicit function using Eq. (3), velocity function values are required at all grid points. A modified fast 
marching method suitable for unstructured grids is used to extend the velocity function values to all grid points [23]. 
This approach also has the benefit of maintaining the signed distance function, which is important for the stability of 
the level-set method. 


III. Aeroelastic Modeling 

This section presents brief details on the modeling strategy used for static and dynamic aeroelastic analyses. The 
subject of the aeroelastic optimization in this paper is an aircraft wing in subsonic flow. 

A. Static Aeroelastic Analysis 

The methods described in this section compute the trimmed static aeroelastic equilibrium of the wing. The 
deformed state can then be used to compute buckling load factor eigenvalues for the wing skins. Stress and stiffness 
quantities measured within this trimmed state can also be considered during the optimization, but this is not done 
here. There are three main numerical tools used to compute the static aeroelastic equilibrium of the wing. FEA is 
used to compute the structural equilibrium: 


K(Q).u = Af.f,(Q)+f4c^) (5) 

where K is the structural stiffness matrix, u the displacement vector, N is the load factor (1 for steady level flight), ft 
is the self-weight body force load vector and fa is the aerodynamic load vector, which is dependent on the pressure 
distribution, Cp , computed using the Doublet Lattice Method (DEM) [26, 27]: 

D-c^ = Wj.-a(u,Q)-z + w(u) (6) 

where D is the non-dimensional aerodynamic influence coefficient (AIC) matrix, Wc is the slope of the aerofoil 
centerline (from built-in camber and twist), z is a column vector of ones, w(u) is the additional downwash dependent 
on the deformed shape of the wing and a is the angle of attack at the wing root computed to meet the trim condition 
that lift equals weight: 


L[c^) = q-a^-c^{a) = N-W{n) (7) 

where L is the total lift from one wing, W is half the weight of the aircraft, q is the dynamic pressure and a is a 
vector of box areas in the DEM discretization. 

The Finite Plate Spline (FPS) method [28] is used to couple the three equilibrium equations, Eqs. (5-7). This 
method communicates loading and displacement information between the FE mesh and DEM discretization through 
another FE mesh composed of plate elements. The communication process can be written using linear equations: 

f,{c^) = q-TpC^ ( 8 ) 

w(u) = T„-u (9) 
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where T/is the load transfer matrix and Tu is the displacement transfer matrix. The data transfer matrices are derived 
by applying a series of displacement constraints at points in the FPS plate mesh that coincide with DLM box centers 
and a sub-set of the FE mesh nodes. 

Combining equations (6-9) the aerodynamic load in Eq. (5) can be written as a function of the displacement 
vector and structural domain: 


f„(u,Q) = f,-a(u,Q)-f„+Q-u 


( 10 ) 


where the angle of attack for trim is written as: 

a(u,Q) = + b • u - ^ • ff(Q)) /L„ 


( 11 ) 


the terms S, C , fa , Q, b, Lc and La are all constants, which are listed below in Eq. (12). Full details on the 
derivation of Eqs. (10-12) are in [23]. The displacement vector for static equilibrium of the aeroelastic wing can be 
found by solving a single linear equation, by substituting Eqs. (10, 11) into Eq. (5). 



(12a) 

= ^ S D ‘ z 

(12b) 

Q = ^ S D ‘ T 

(12c) 

b^=q D ' T 

(12d) 

4 = ^-a^-D-‘-w, 

(12e) 

N 

b 

II 

(12f) 


B. Linear Buckling Analysis 

A linear buckling eigenvalue analysis is used to compute the buckling load factors [29] using the displacement 
vector for the trimmed aeroelastic condition computed in the previous section: 

[K(Q,t) + A-K^(u)]-8 = 0 (13) 

where is the stress stiffness matrix, ^ is a eigenvalue corresponding to a buckling load factor, 5 is an eigenvector 
corresponding to the buckling mode shape, u is displacement vector from the static aeroelastic equilibrium Eq. (5). 

In this work we only consider buckling in the wing skins, modelled using shell elements and buckling is not 
considered in the internal structure described by the level-set function. Therefore, the stress stiffness matrix is only 
non-zero for degrees of freedom associated with the wing skins. This approximation was compared with the case 
where is computed for all degrees of freedom and the resulting buckling load factors were within 0.5%, thus 
validating the approach. The shell thickness values can also be design variables and optimized simultaneously with 
the internal structure using the SEP level-set method. However, in this paper we use fixed shell thickness values and 
only optimize the internal structure described by the level-set function. 

A buckling load factor can be written as a function of the structural domain and aeroelastic displacement vector 
by pre-multiplying Eq. (13) by 8^ and rearranging: 


A(Q,u(Q)) = - 


6^K(Q)-6 

8^-Kj(Q,u(Q))-8 


( 14 ) 
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C. Dynamic Aeroelastic Analysis 

The dynamic aeroelastic phenomenon considered in this paper is flutter. The p-k method is used to compute the 
flutter point. The p-k flutter equation used here has an aerodynamic damping term derived from the imaginary part 
of the unsteady aerodynamic matrix [30]. Therefore, the p-k equation is a quadratic eigenvalue problem, which can 
be written in first order form as: 


F[U,q,MX^) 


I 0 

0 (f)'M(Q) 


0 I 


V 

P'\ 


= 0 


(15) 


where p is an eigenvalue, M is the structural mass matrix, Ar and Ai are the real and imaginary parts of an unsteady 
aerodynamic matrix, respectively, where the aerodynamic matrix is non-linearly dependent on the reduced 
frequency k and Mach number M, v is an eigenvector, U is the air speed and Z? is a reference length, taken here as 
half the root chord. The aerodynamic matrices are obtained using the DLM and FPS method: 

A(M,A:) = T,-D-'(M,A:)-(T„ + /-fT,) (16) 

where Tv is a matrix computed by the FPS method that determines the downwash due to the wing motion normal to 
the planform from the structural velocity vector. The AIC matrix used in Eq. (16) is evaluated for A: > 0 and is 
complex, whereas the AIC matrix used in the static aeroelastic analysis, Eq. (6) is evaluated for A: = 0 and is real. 

A key aspect of the p-k method is that the input reduced frequency k must equal the imaginary part of the 
resulting eigenvalue: Im(p) = A:, and the air speed and Mach number must also correspond for a given altitude. When 
these two criteria are satisfied then a matched point solution to Eq. (15) has been found. The real part of the 
eigenvalue for a matched point approximates the required damping to make that particular mode stable. If Re(p) < 0, 
then that mode is stable. If Re(p) > 0, then some damping is required to stabilize the mode. If all of the modes are 
stable, than the complete system (wing) is stable. The flutter point is often defined as the matched point solution 
with the lowest air speed or dynamic pressure that produces a single eigenvalue with Re(p) = 0, with all other modes 
stable. 

The matched flutter design point can be written as two equations that are dependent on the solution of Eq. (15): 


mp) K 

lm.{p)-k J 


(17) 


The flutter point is usually found by increasing dynamic pressure, or air speed, until both conditions in Eq. (17) are 
satisfied by a solution to Eq. (15). The dependence of the aerodynamic matrices on Mach number and input reduced 
frequency, k, is non-linear and contained within the DLM numerical procedure. Both of these non-linear effects 
should be taken into account when searching for a matched point solution to Eq. (15). This usually requires the 
computation and storage of a sufficient number of aerodynamic matrices at various k and M values to model the 
entire non-linear space via interpolation. This computational burden can be reduced if the Mach number is fixed 
during the analysis. If the maximum air speed is sufficiently small compared to the speed of sound at the design 
altitude, then compressibility effects can be ignored and the Mach number fixed at zero. If this is not possible, then 
an alternative is to fix both the air speed and Mach number in Eq. (15) and define the flutter point in terms of 
dynamic pressure or air density. This approach is equivalent to defining the flutter point as the minimum altitude at a 
given speed where: Re(/7) = 0. This approach is still an approximation, as the speed of sound is generally not 
constant with altitude, although the dependence is relatively weak. In this work, the Mach number and air speed are 
fixed and the flutter point is defined as the minimum dynamic pressure that produces an eigenvalue solution to Eq. 
(15) that satisfies Eq. (17). 

The flutter point is first approximately found by incrementally increasing q from zero until the real part of an 
eigenvalue becomes positive. The exact matched point is then found using Newton’s method to satisfy the residual, 
Eq. (17) [7]. Matched points in the k space at a specific q value are found using a non-iterative frequency sweeping 
technique [31]. For a set dynamic pressure, k is incrementally increased and the corresponding eigenvalue p is 
monitored until a crossing with the unit-slope line (in the Im(p) - k space) is noted. The matched point solution for 
k = Im(p) is then obtained by linear interpolation. The sweeping approach is more efficient than frequency iteration 
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methods, as it requires fewer solves of Eq. (15). It is also more robust, as the frequency iteration approach is not 
guaranteed to converge to a matched point. A critical aspect of the frequency sweep approach is the tracking of 
eigenvalues as k is increased (for a given q) and also as q is increased. Mode tracking is performed using an 
eigenvector orthogonality correlation method [32]. The correlation matrix is formed from scalar products between 
the current eigenvectors and the converged eigenvectors from the previous step. The values in the matrix are then 
used to re-correlate the modes. 

The order of the dynamic aeroelastic system in Eq. (15) for a full 3D wing box is large and the eigenvalue 
analysis is computationally expensive. A reduced order method using the first 15 free vibration modes of the 
structure is used to decrease the degrees of freedom in Eq. (15) to a small number of generalized degrees of freedom 
[29]. For example, the reduced mass matrix is: 


(18) 

where Mr is the reduced mass matrix and ^ is a matrix where the i^^ column contains the i^^ eigenvector 
corresponding to the i^^ free vibration mode of the structure. Reduced stiffness and aerodynamic matrices are 
obtained in the same way. The reduced order matrices are used in place of the full order matrices in Eq. (15) to more 
efficiently compute the flutter point. 


IV. Aeroelastic Optimization 

This section describes several static and dynamic aeroelastic optimization problems and sensitivities are derived 
for the aeroelastic objectives and constraints. 

A. Buckling Problem 

The first problem is to minimize the mass of the wing subject to a constraint that buckling does not occur under 
static aeroelastic loading: 


Minimize: Mass 
Subject to: >1.0 


(19) 


where Anin is the minimum buckling load factor. Since the solution to the buckling eigenvalue problem clearly 
depends on the static aeroelastic trim state, Eq. (13), the adjoint gradient method must reflect this dependence. The 
total derivative of the buckling load factor, Eq. (14), with respect to the structural domain is: 


dX dX T ^(^) r ^.(^) 

— = — -p ^-u + p • ^ + ^ 

dQ tTl ^ ^ tdn _ 


where: 


dX J.X 5K(Q) 5 . 

— = -o • o 

5f _ da(Q) ^ _ N dW(Q) ^ 
5Q 5Q " z/ 5Q ' 

and p is the solution to the adjoint equation: 



'^7 

T 


p = 


p + 

-Al. 


( 20 ) 


( 21 ) 

( 22 ) 


( 23 ) 


where: 
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( 24 ) 


du 

J _ f„®b 

5u 5u “ 

The constant aerodynamic terms in Eq. (25) are defined in Eq (12). Note in Eqs. (21, 24) the eigenvector 8 has been 
normalized to make the denominator in Eq. (14) equal to one. The derivative of the stress stiffness matrix with 
respect to displacement in Eq. (24) is obtained by differentiating the shell element matrices used to form K^. 

There are many solutions to the linear buckling eigenvalue problem, Eq. (13), which represent many buckling 
modes. Many of the load factors associated with these modes may become critical during optimization. The 
constraint in Eq. (19) is for the minimum load factor to be greater than one. The switching of critical buckling 
modes during optimization presents a potentially discontinuous design space, which may present convergent 
difficulties in gradient-based optimization. This can be remedied by specifying separate constraints for a sufficient 
number of modes associated with the most critical buckling load factors. However, this results in a large number of 
constraints, which significantly increases the computational cost of the optimization. A common approach to deal 
with a large number of constraints is to use an aggregation technique to convert many constraints into a single 
constraint. The Kreisselmeier-Steinhauser (K-S) function [33] is used here to aggregate a number of the most critical 
buckling load factors: 




Jexp(r-(A^„-A,)) 


(26) 


where r is the aggregation parameter and n the number of buckling modes used in the aggregation. Modes with the 
lowest positive buckling factors are included in the aggregation. The K-S function provides a conservative estimate 
of the original constraint. The buckling optimization problem with constraint aggregation then becomes: 


Minimize: Mass 
Subject to: KS{X) < 0 


(27) 


The derivative of the K-S function is composed of derivatives for each constraint function and only requires a single 
adjoint solution [33]. 

B. Flutter Maximization Problem 

The optimization problem is to maximize the flutter dynamic pressure subject to a mass constraint: 


Maximize: 

Subject to: Mass < Mass * 


(28) 


where is the matched flutter dynamic pressure andMa^^* is the mass upper limit. The derivative of q^ is obtained 
by taking total derivative of the residual equation that defines the matched flutter point, Eq. (17), at the matched 
point, setting it to zero (as the residual will still be zero at the next iteration) and rearranging: 

I dq^/dO. 

I dk^/dQ. 

where A:* is the reduced frequency at the matched point and the Jacobian matrix in Eq. (29) is obtained at the end of 
the Newton’s iterations used to locate the exact matched point [7]. The derivative of the eigenvalue p is obtained 
using derivatives of the p-k flutter eigenvalue problem, Eq. (15): 


^Re(p) 

dq* 

^Re(p) 

dk* 

-1 

< 

Re(5/>/5Q) 

dlm{p) 

dq* 

dlm{p) 1 
dk* ^ 


lm(dp/dQ!^ 


(29) 
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( 30 ) 


The derivatives of Eq. (15) are: 




dp 


dp 

dF jdF 
dQ./ dp 


0 

0 

n 


L, r 

I 0 

1 R. 1 


1 

% 

O 

1 

|r. J 



(31) 


(32) 


where L and R are the left and right eigenvectors at the matched flutter point and superscript H is the conjugate 
transpose. To compute the shape derivative of the flutter dynamic pressure we require the derivative of the reduced 
stiffness and mass matrices. To achieve this we start by multiplying out each entry in the reduced matrices: 


K=WK^- 


^1,1 

^1,2 




mi2 • 


^2,1 

^2,2 

*• ^2,n 

, = 

'”2,1 

^2,2 * 




Cn,n 






(33) 


where n is the number of structural eigenvectors used to form the reduced order model. The matrix entries are 
approximated using a continuous weak form: 




LEe{Vi)-4Wj)d^ , m,j^\ p-it/py/jdQ. 


(34) 


where E is the material property tensor, s the strain tensor, p the material density and y/i is the i^^ natural frequency 
eigenvector. The shape derivatives of the functions in Eq. (34) are well known [34] and can be used with Eqs. (29- 
34) to obtain the shape derivative of the flutter dynamic pressure. Fixed mode derivatives are used that ignore the 
derivative of the eigenvectors with respect to Q when computing the shape derivative, although the mode shapes are 
updated from one design iterate to the next. This is a common approximation, as the computation of eigenvector 
derivatives is expensive [35]. 

C. Flutter Constraint Problem 

The problem is to minimize the wing mass subject to a constraint that the flutter dynamic pressure is greater than 
a limit value: 


Minimize: Mass 
Subject to: ^* > 


(35) 


where qum is the limit on the flutter point dynamic pressure. 

The flutter point is often a non-linear function of the design variables and it can also be discontinuous, due to a 
switching of the critical flutter mode. To avoid problems caused by a discontinuous constraint function, the flutter 
constraint is reformulated as a parametric constraint [36]. The optimization problem is then: 

Minimize: Mass 

Subject to: G^{q) < 0 , <q<q^^ 
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where qmin and qmax are limits on the dynamic pressure range of interest and G/is a continuous function defined as: 



if 

& Rq(p)>0 

P,ef 

if 

& Re(/?)<0 


1 if 

& Re(/?)>0 


if 

& Re(/?)<0 


(37) 


where pref is a reference value, defined here as the reference length divided by the air speed: pref = b I U. Note that 
the third component of Eq. (37) (for q < qum & Re(p) > 0) has been added to the originally proposed formulation in 
[36] to ensure the constraint function is continuous for all possible q and Re(p) values. 

The p-k method and Newton’s method are used to find local maxima of G/ in a similar way that the matched 
flutter point is found, except the Re(/7) = 0 condition is not required. Each local maxima indicates a point that could 
become the critical flutter mode and therefore must be accounted for in the optimization problem. Thus, although the 
problem in Eq. (36) is stated with one constraint, G/ can have several local maxima, leading to several constraint 
functions that must be handled by the optimizer. This is handled by the SEP level-set method by formulating the 
SEP sub-problem used to obtain the velocity function with as many constraint functions as there are local maxima. 
However, if a local maximum occurs far from the Re(p) = 0 line, it is not critical and thus ignored in the SEP sub- 
problem. 

The shape derivatives of local maxima of the continuous constraint function require the total derivative of Re(p): 

JRe(/>) _ 5Re(/>) dRQ(p) dq dRQ(p) dk 
dQ ~ 5Q ^^q ^ ^ 


For a small change in the design the continuous constraint function can be approximated about the current q value 
and the dq/dCl term in Eq. (38) can be assumed zero [36]. The first term on the right of Eq. (38) is obtained from 
Eqs. (30-34). Taking the total derivative of Eq. (15) with respect to k, setting it to zero and rearranging gives: 


dp _ dF j dF 
dk dk / dp 


(39) 


where the denominator is obtained from Eq. (32) and the numerator is obtained by taking the partial derivative of 
Eq. (15), pre-multiplied by the complex conjugate of the left eigenvector at the local maximum. This involves 
computing the derivatives of the real and imaginary parts of the reduced aerodynamic matrices with respect to k. 
These are obtained using linear interpolation of the pre-computed set of aerodynamic matrices at specific k values 
that are used during the p-k method. The final component required for Eq. (38) is obtained by taking the total 
derivative of the residual equation for a matched point: Im(/7) - k=0, setting it to zero and rearranging: 

^ ^ aim(/>) jL_ aim(j9) "| 

da XI F dk ) ^ 


where the numerator is obtained using Eqs. (30-34) and the denominator is obtained at the end of the Newton’s 
iterations used to find the exact matched point. 


V. Examples 

The SEP level-set method is used to solve several aeroelastic optimization problems for a realistic wing. The 
geometry is derived from the Common Research Model (CRM) wing [37]. The CRM outer mould line (OML) 
represents the cruise shape of the wing, under aerodynamic loading. For aeroelastic analysis the structural mesh 
should represent the unloaded shape of the wing, commonly referred to as the jig shape. Therefore, a representative 
jig shape OML is used in this study. The CRM is a transonic transport configuration with an aspect ratio of 9, a taper 
ratio of 0.275, a sweep angle of 35° and a cruise Mach number of 0.85. The wing box is located between 10% and 
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70% of the local airfoil chord. The wing box structure is made from aluminum, with a Young’s modulus of 68.95 
GPa, Poisson’s ratio of 0.3 and density of 2800 kg/m^. 

The internal structure of the wing box is represented by the level-set implicit function, Eq. (1), discretized on a 
3D grid fitted to the wing box shape. This part of the structure is optimized by iteratively solving Eq. (3), where the 
velocity function is obtained from shape sensitivities and an SEP sub-problem. The upper and lower skins and 
leading and trailing edge spars are modeled using shell elements that are connected to the solid elements. Lumped 
mass elements are used to model the effect of leading and trailing edge devices and an engine. The mass elements 
are attached to the wing box structure using an interpolation method that distributes internal forces such that no 
additional stiffness is introduced [38]. The solid elements are isoparametric 8 node bricks with incompatible modes 
to improve bending performance [29]. The shell elements are 3 node facet shells, each composed of a discrete 
Kirchoff plate element and a linear strain triangle, with mid-side degree of freedoms replaced by drilling degrees of 
freedom at the corner nodes [29]. 

The solid element mesh for the CRM wing box is shown in Fig. 1. The full FE model has 154,560 brick 
elements, 59,248 shell elements, 18 mass elements to represent the leading and trailing edge devices, and an 
assembly of 13 mass elements to model the engine, shown in Fig. 2. The 12 outer engine mass elements are 
connected to a central element by rigid links. The central element is then connected to the leading and trailing edge 
spars using the interpolation method. The model has a total of 170,252 nodes and 597,954 degrees of freedom. All 
degrees of freedom for nodes at the wing root are fixed. The DEM mesh is composed of 20 x IQO boxes in the 
chord- wise and span-wise directions, respectively. A subset of 1665 nodes on the upper surface and 1665 
corresponding nodes on the lower surface are used with the EPS method to compute the transfer matrices. The 
complete static aeroelastic model and p-k method, described in Section III, were validated against results obtained 
from MSC NASTRAN. 

The mesh used to discretize the implicit level-set function is reasonably coarse. Therefore, thin shell-like 
structural members, typical of a conventional semi-monocoque configuration, are not easily represented inside the 
wing box. Therefore, the structures obtained using the mesh in Fig. 1 are likely to be much heavier and stiffer than 
real wing boxes. To obtain more realistic results, the modulus and density of the material used for the internal 
structure are a fraction of the real values, which is specified in each example below. 



Figure 1. Solid mesh of the CRM wing box. 
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Figure 2. CRM mass element attachments. 


A. Buckling Constraint Problem 

The buckling constraint problem, Eq. (27), is solved using 75 of the lowest positive buckling load factors in the 
K-S function. The material properties used for the internal structure are 20% of the real aluminum values. The shell 
thickness values vary linearly along the span from 15 mm at the root to 4.5 mm at the tip. A single aeroelastic load 
case is considered that corresponds to a 2.5g pull-up maneuver (N= 2.5 in Eqs. (5, 7)). The dynamic pressure is 8.8 
kPa and the weight of the aircraft, excluding the wings, is 500 kN. The initial structure is a series of ribs with thicker 
ribs at the tip, root and where the engine is attached, as shown in Fig. 3. The solution is shown in Fig. 4, where the 
initial ribs have been reduced to chord-wise stiffeners located in the same place as the initial ribs. This suggests that 
the solution may be dependent on the initial design. The convergence history is shown in Fig. 5, which shows an 
initial increase in the K-S buckling function, followed by the constraint becoming feasible, oscillating about the 
limit value of zero for about 30 iterations, before reaching a converged solution after 68 iterations. The minimum 
buckling load factor for the solution is 1.07. This slight conservatism is introduced by the K-S constraint 
aggregation, Eq. (26). This example shows that the SEP level-set method can obtain feasible local optimum solution 
for a buckling constrained problem considering aeroelastic trim loads. 



Figure 3. Buckling example initial internal structure: a) 3D view, b) top view. 
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Figure 4. Buckling constraint problem solution. 
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Figure 5. Buckling constraint problem convergence history. 

B. Flutter Maximization Problem 

The flutter dynamic pressure maximization problem, Eq. (28), is solved with a fixed air speed of 250 m/s and 
constraint of 11,161 kg on the total wing mass, which corresponds to a 50% volume fraction for the internal 
structure. For this example 1% of the real aluminum values are used for the internal structure and shell thicknesses 
are 2 mm. The initial internal structure is shown in Fig. 6 and the solution is shown in Fig. 7. The solution has 
several interesting features. Firstly, the outboard section leading edge region is reinforced with material, whereas the 
trailing edge region is mostly void. This moves the center of gravity forward, which is a mass balancing effect 
known to be favorable for flutter [39]. There is some lattice type structure near the root, which is evidence of an 
aeroelastic tailoring effect, although the exact role of the structure is unclear. 

The convergence history is shown in Fig. 8. Switching of the critical flutter mode causes the discontinuities at 
iterations 24-25 and 42-44. The damping curves for iterations 42-44 show that the critical mode switches from mode 
2 to mode 13 at iteration 43, then back to mode 2 at iteration 44, Fig. 9. The critical flutter displacement mode 
shapes for iterations 42-44 show that the critical mode for iterations 42 and 44 are global modes involving motion of 
the entire wing. Fig. 10. However, the critical mode for iteration 43 is a local mode contained within a single panel 
near the root. Fig 10b. The local mode in this model is spurious because the DFM panel method assumes that the 
upper and lower surfaces move in tandem, which is clearly not the case for the local mode. One method to avoid the 
appearance of spurious local panel modes is to place spline points on the FE mesh at “hard point” that coincide with 
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ribs and spars and not in the middle of a panel [40]. However, for topology optimization of a wing-box there are no 
set hard points as the positions and number of internal member can change during optimization. An alternative is to 
handle the mode switching using the continuous formulation employed for the flutter constraint problem, Eq. (37), 
which will prevent local modes from becoming critical. 



Figure 6. Flutter example initial internal structure: a) 3D view, b) top view. 



Figure 7. Solution for flutter maximization problem: a) 3D contour of internal structure, b) sections 

though the contour (not ribs). 
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Figure 8. Convergence history for flutter maximization problem. 
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Figure 9. Flutter maximization problem damping curves: a) iteration 42, b) iteration 43, c) iteration 44. 



Figure 10. Critical flutter mode shapes: a) iteration 42, b) iteration 43, c) iteration 44. 

C. Flutter Constraint Problem 

The dynamic pressure limit in the continuous constraint function, Eq. (37), is specified as 38.2 kPa, which equates 
to approximately sea level air density using a fixed air speed of 250 m/s. For this example 5% of the real aluminum 
values are used for the internal structure, the skin thickness is 2.5 mm and the spar thickness is 15 mm. The initial 
design is the same as the previous example, shown in Fig. 6. 

The solution is shown in Fig. 11 and the convergence history is shown in Fig 12. The optimizer removes almost 
all the internal structure. Reinforcement is seen in areas near the leading edge, break and engine attachment points. 
The convergence history shows that the final solution at iteration 136 is feasible as the maximum value of G/ is 
negative. The design does not stay feasible throughout the optimization and oscillations start to occur after iteration 
35. The oscillations are mainly caused by a hump mode becoming critical from small changes in the design. This is 
highlighted by the damping curves for iteration 131 and the final design, shown in Fig. 13, where modes IB, IT and 
IE correspond to the P* bending, P* torsion and P^ edgewise bending modes, respectively. This mode labeling is 
based on the dominant component of the mode, as the natural modes shapes of the final design represent 
combinations of modes. For example Fig. 14 shows that the critical flutter mode shapes (P* torsion and P* edgewise 
bending) also have a significant bending component. The continuous constraint formulation, Eq. (37) should prevent 
the discontinuous switching of modes seen in this example. It is thought that inaccuracies in the sensitivities lead to 
an underestimation of how the hump mode local maximum changes as the design is updated, resulting in the 
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observed mode switch, although this requires further investigation. Inaccuracies may come from the fixed mode 
derivative assumption. 



Figure 11. Solution for flutter constraint problem: a) 3D contour of internal structure, b) sections though 

the contour (not ribs). 



Figure 12. Convergence history for flutter constraint problem. 



Figure 13. Flutter constraint problem damping curves: a) iteration 131, b) iteration 136 (final). 
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Figure 14. Mode shapes for final design: a) 1*^ torsion, b) 1*^ edgewise bending. 

VI. Conclusions 

This paper presents a strategy for using level-set topology optimization to solve problems that consider static 
aeroelastic loads and dynamic aeroelastic stability. The SLP level-set method is used to obtain a velocity function 
that improves the objective, whilst maintaining constraint feasibility. The proposed method is used to optimize a 
representative CRM wing box. The internal structure of the wing box is modeled using an implicit level-set function 
discretized on a fitted 3D grid, whereas the skins and leading and trailing edge spars are modeled using shell 
elements. This modeling approach allows skin buckling to be optimized effectively. Lumped mass elements are used 
to model leading and trailing edge devices and an engine. The proposed method is demonstrated to solve three 
problems: mass minimization with a linear buckling constraint, mass minimization with a flutter constraint and 
flutter maximization with mass constraint. A K-S function is used to aggregate many buckling constraints into a 
single constraint function that provides a conservative estimate of the original constraint. A continuous constraint 
formulation is used to manage the discontinuities in the flutter constraint, such as switching of the critical mode. The 
strategies used were effective at obtaining feasible optimal designs, although it appears that solutions may be 
dependent on the initial design. 
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